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Abstract 

The phase diagram and critical behaviour of a simple toy model for DNA zipping/unzipping is 
examined in the framework of the Bethe approximation. The effects of solvent quality are included, 
and found to lead to a variety of different thermodynamic behaviours. 
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I. INTRODUCTION 



Self-avoiding walk models and their variants have been used for decades to gain insight 
into the physics of real polymer systems, with remarkable success 111]. Such models are now 

nnnnn 

again of interest as toy models of biopolymers such as DNA and RNA molecules |2J, y|, |4], |5|, la] ■ 
The effect of the base pairing is modelled by allowing the walk to visit the lattice bonds 
twice. In the context of the work presented here, the difference between DNA-like models 
and RNA-like models are the allowed walk configurations; DNA models allow at most four 
links of the walk to meet at a site (see figured]), whilst RNA allow up to eight. 

Recently several models of this type have been presented, either as models for biopoly- 
mers, or simply as interesting variants of restricted random walk models, in which different 
weightings are given to multiply visited sites or bonds 3, Isl, What is missing from all 
these models are interactions representing the interactions with solvent molecules. In this 
article we propose a first step at filling this void. 

In this article, we introduce a variant of the lattice two-tolerant self-avoiding walk 
model 3, 10] which mimics the zipping/unzipping of a DNA model. In the model in- 
troduced here we also include solvent effects by including an attractive interaction energy 
between non-consecutively visited nearest-neighbour sites. In this article we present an ex- 
tended mean-field type calculation of the phase diagram (the Bethe approximation). The 
phase diagram is found to be unexpectedly rich. In many DNA models there is a further 
restriction in that bases are only allowed to pair up if they are the same distance from one 
end along each of the two chains. This restriction is relaxed in the model presented here. 



II. MODEL 



The model studied in this article consists of a non-crossing random walk on a square 
lattice limited to visit each bond and each site at most twice. The model is chosen to 
model the zipping/unzipping of a DNA molecule, thus the allowed configurations of the 
two-tolerant walk are further restricted such that a site may only be visited twice if one of 
the adjoining bonds is doubly visited. The allowed configurations are shown in figure HJ 

Each segment of the walk represents a monomer (base) whereas doubly-visited bonds 
represent a coarse-grained description of paired bases. The difference in affinity of the DNA 
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FIG. 1: Configurations at a lattice site on the square lattice. Polymer bonds are represented by solid 
lines while the solvent interactions are represented by dashed lines. The empty-site configuration 
is not shown here. 

molecule with itself and with the solvent may be modelled by a solvent-mediated attractive 
interaction between non-consecutive visited nearest-neighbour sites on the lattice. Solvent- 
mediated interactions carry an attractive energy — es and doubly-visited bonds yield an 
attractive energy —e. 

The thermodynamic behaviour may be investigated by introducing the grand-canonical 
partition function, Z, from which many of the relevant thermodynamic quantities may be 
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calculated. The grand-canonical partition function is given by: 

z=J2k n t Ni l N2 (1) 

walks 

where Nj is the number of solvent-mediated interactions, N 2 is the number of doubly visited 
bonds and r = exp (— Pes), 7 = exp (—/3e) and (3 = 1/kT. The fugacity, which controls the 
average length of the walk, is denoted by K, and N is the total length of the walk. The two 
are related through: 

W (2) 
The average length increases as K is increased. 



III. THE BETHE APPROXIMATION 

In this section we briefly describe the Bethe approximation. For a good discussion of 



the Bethe approximaton see Ref . . The model of interest is studied on the infinite Bethe 
lattice chosen to have the correct local geometry. The lattice chosen for the square lattice is 
shown in figure [2] The Bethe lattice is a hierarchical lattice built recursively from a central 
bond by adding k new bonds to each extremity. To each dangling bond we add k more 
bonds, and so on, such that no loops are formed. Due to the hierarchical nature of the 
lattice, it is possible to build up expressions for the partition function recursively. To see 
this, it is convenient to consider the lattice as being divided into two branches, left and 
right for the example shown in figure [2] We may introduce the partial partition functions 
W\ and W T a for the left and right-hand branch, respectively. These partition functions are 
conditional upon the state a of the central bond. In our model there are four possible states: 

1. empty (state 0), 

2. occupied with a link of the walk (state 1), 

3. occupied with a solvent-mediated nearest-neighbour interaction (state S), 

4. occupied with a doubly- visited bond (state 2). 

By symmetry, the left and right branches will have the same partial partition functions, 
and so the l,r designation will be dropped. Each branch may be sub-divided into k sub- 
branches, such that the W a may be expressed in terms of the partial partition functions of 



the sub-branches. This procedure may be continued until the boundary bonds are reached. 
In order to do this explicitly, it is convenient to introduce the notion of the 'generation' of a 
link, n, which is simply the distance of the link from the boundary. As a concrete example, 
consider the calculation of the partial partition function conditional on the central 

bond being occupied by a link of the walk. We must consider all the configurations on the 
bonds of the generation (n — 1), of which there are three for the 2 dimensional square lattice 
example shown in figure [21 which are compatible with the occupied central bond. Clearly 
there must be a bond leaving in one of the three directions, the other two bonds may be 
empty or occupied by a solvent-mediated interaction. The weight is simply the sum of 
the Boltzmann weights corresponding to all these configurations, multiplied by the weight 
for adding the central link. To avoid the divergence of the partial partition functions it is 
convenient to introduce normalised partition functions = W^/q n 4\ with q n chosen 
such that: 

= (3) 



FIG. 2: The Bethe lattice representation of the two-dimensional lattice. The dotted box shows the 
central bond, exhibiting the desired square-lattice geometry. 
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This leads to recursion relations for the (normalised) partial partition functions: 

# = ^E^,wn< -l) > (4) 

{a,} 1=1 

where {<7j} is the set of states of the k links forming generation n — 1, A CT is the Boltzmann 
weight of the bond added at generation n, and the factor C^o-j = 1 if the choice of the 
states {<7j} is compatible with the central state a, and zero otherwise. 

It is known that there is no phase transition on the infinite Bethe lattice, since the number 
of boundary sites grows too rapidly. However the recursion relations may be used in the 
centre of the lattice as self-consistency equations for the two point mean-field theory for the 
corresponding square lattice. In this case, we assume we have translational invariance, and 
drop the generational superscripts. The equilibrium states are then given by solutions of 
the following set of recursion relations: 



w = - [u>o + 3w\ (w + w s + w 2 ) + 3u>2 (w + w s )j (5) 
K 

w 1 = 3— w 1 (w + w s ) (w + w s + 2w 2 ) (6) 

w s = 3 ^ T - — {w\ (w + w s + w 2 ) + w\ (w + ^s)) (7) 

w 2 = 3— — (w + w s ) (w\ + w 2 (w + w s )) (8) 
q = Wq + 3wj (w + ws + w 2 ) + 3wj (w + w s ) 

+3Kwi (w + w s ) (w + w s + 2w 2 ) (9) 
+3fC 2 7 (w + w s ) (w\ + w 2 (w + w s )) 



+3 (r - 1) [w\ (w + w s + w 2 ) + w\ (w + w s )) ■ 

The partial partition functions give the contribution to one branch of the total partition 
function, the total (normalised) partition function conditioned upon the state of the central 
bond is then given by the product of the weights for the left and right branches. Each of the 
partial partition functions includes the Boltzmann weight corresponding to the state of the 
central bond, which is thus counted twice in the full partition function. This double counting 
is corrected by dividing each term by the relevant Boltzmann weight. Summing over all the 
possible states for the central bond gives the total (normalised) partition function, z: 

^ = E^- (io) 
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In the usual way, the probability of finding a given bond in state a is given by the partition 
function conditioned upon this state divided by the total partition function, i.e. 

w 2 

It should be noted that the density p of the walk on the lattice is simply 

p = Pl + 2p 2 (12) 
Another quantity of interest is the fraction of paired segments, given by 

$ = 2 -f (13) 

The grand potential per site may be related to z and q through the relation 

(k - l)\nz - 21ng 
Pf=- — =\nz-lnq, (14) 

for the square lattice (k = 3). For a full derivation of this expression see 4j]. When multiple 
solutions to the recurrence relations exist, the solution with the lowest value of the grand 
potential is the stable equilibrium solution. 

IV. RESULTS 

It is convenient to recast the recursion relations 0— M by setting: 

Wl 

Xt = — , 
w 

w s 

xs = — , 

Wo 

w 2 

x 2 = — . 

Wo 

The different phases and transition lines correspond to different solutions of the recursion 
equations. These equations have a trivial solution x\ = x$ = x 2 = 0, corresponding to 
the finite polymer phase. The density is trivially zero, since the lattice is infinite. In what 
follows we refer to this phase as the O phase. Non-trivial solutions may be found by setting 
different values of the parameters K, t and xs and then solving numerically for 7, x\ and 
x 2 . 

The trivial zero-density phase is separated from a region where the walk fills the lattice 
with a finite density by a surface K = K 00 {r 1 r y) where the average length of the walk first 




FIG. 3: The phase diagram in the K — 7 plane, with r = 1. O denotes the zero density phase, I 
the ordinary dense phase and II the fully paired dense phase. Dashed and solid lines denote first- 
and second-order transitions respectively. The tricritical point on the O-I transition line is given 
by K = 1/3 and 70 = 9/5. The second-order transition between phases O and II is given by the 
line 3-ftT 2 7 = 1, and finishes in a critical end point. The dotted line is the extension of this line and 
is included to guide the eye. 

diverges. In the finite-density region we find two dense phases : a phase (I) in which a finite 
fraction of the segments of the walk are paired, i.e. it is characterised by < p < 1 and 
< $ < 1 ; and a fully paired dense phase (II) in which every segment is paired i.e. it is 
characterized by < p < 1 and $ = 1. 

When 7 = (giving rr 2 = 0) the pairing of segments is forbidden, and the model cor- 
responds to the standard interacting self-avoiding walk model. For this model, also called 

n 

the 0-point modelpj, phases O and I are separated by a critical transition line for small 
enough r and a first order line for large enough r. These two behaviours are separated by 
the tricritical O-point. This tricritical point extends to a line of tricritical points as 7 is 
allowed to increase, separating a critical region in the SAW universality class from the first 
order region, where the walk fills the lattice with a finite density. It is possible to determine 
analytically the region occupied by the SAW phase, and so the equation of the tricritical 
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line, as follows: from phase I (x± 7^ 0) we take the limit x$ — > and £2 — > in order to 
approach the boundary with phase O. Since the transition is continuous, this limit leads to 
x\ — > and gives = l/(2d — 1) = 1/3. The tricritical point is obtained by allowing 
the parameters x\, xs and X2 to be vanishingly small but non-zero along the transition line 
K = 1/3. We obtain the tricritical condition : 

9(3 - 2r) 



7e (r) 



(15) 



4(r- l) + 5(3-2r) 

It is important to note, that there is no guarantee that the extended line is in the same 
universality class as the O-point, but we will nevertheless use the B subscript to differentiate 
it from other special points in the phase diagram. We shall return to this point later. If the 
pairing interaction is taken to be attractive, then 7 > 1, this leads to the result that the 
transition line O-I becomes fully first order for r > 4/3. 

T 




FIG. 4: The two order parameters, the density p and the fraction of paired segments $ as a function 
of 7 for t = 1 and K = 0.278. Solid and dashed lines denote their values in phases I and II. 

The phase diagram for r = 1 is shown in figure [3l clearly showing the location of the 
tricritical 0-like point, found to be K = 1/3 and 70 = 9/5. The phase boundary of phase 
I is smooth, with the boundary between phases I and II being partly first order and partly 
second order. This gives rise to another tricritical point at 7 tc w 4.7. The transition line 
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FIG. 5: The phase diagram in the K — 7 plane, with 4/3 < r = 1.4 < 3/2. Dashed and solid 
lines denote first- and second-order transitions respectively. The point marked cep is a critical end 
point, whilst tc indicates a tricritical point. The dotted line corresponds to the line 3K 2, y = 1. 

between phases O and II is found to be continuous, and the equation for this line may be 
determined exactly: from phase II (x\ = and £2 7^ 0) we take the limit xs — > to approach 
the boundary with phase O. Using the fact that the transition is continuous, we determine 
that £2 — > leads to the condition 3K 2/ y = 1. This line terminates in a critical end point at 
7ce P « 3.94. 

For t = 1, i.e. in absence of solvent effects, the phase diagram is similar to the phase 



diagram obtained by Pretti for an RNA-like model with non-zero stacking energy [4 1 , except 
that in Pretti's case the equivalent to the O-II line was found to be first order. The stacking 
energy had the effect of favouring the absence of multiple double bonds meeting at a site, 
and so made the model more like the model presented here. 

Before examining other values of r, it is of interest to discuss the different phases shown in 
figure El and the relevant order parameters. Two parameters are of interest for differentiating 
the different phases: the density of occupied lattice bonds p, and the fraction of paired 
segments $. These two quantities are plotted in figure H] for K = 0.278, so as to pass 
through each phase. In phase O p = trivially since the walk is finite on an infinite lattice. 
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FIG. 6: The phase diagram in the K — 7 plane, with 3/2 < r = 1.6 < 1.75. Dashed and solid 
lines denote first- and second-order transitions respectively. The point marked tp is a triple point, 
whilst cep denotes a critical end point. The dotted line corresponds to the line 3iT 2 7 = 1. 

In the framework of the Bethe approximation it is not possible to calculate <3> in this phase 
(p = 0,^2 = 0). The value of K is chosen to cross the transition from phase O to phase 
II under the critical end point. We clearly see that p increases smoothly from 0, indicative 
of a critical transition, and that $ = 1 indicating that phase II corresponds to a saturated 
doubly occupied phase. Again we remind the reader that the apparent jump in $ is simply 
related to the fact that it is not possible to calculate $ in phase O (where pi — P2 — 0). The 
value of K is chosen to cut back into phase I between the critical end point and the second 
tricritical point. We can clearly see in figure H] that $ drops disco ntinuously from 1 and 
there is a small jump in p. The transition from phase I to phase II is again continuous, and 
the interest of $ as order parameter is clearly seen; there is no evidence of a phase transition 
in p, which is not a good order parameter for phase II. 

As r is increased, the phase diagram changes. The first change is that the B-like tricritical 
transition is pushed out of the domain of interest (e > or 7 > 1). This occurs, as previously 
stated, for r = 4/3. This case is shown in figure [5j The calculation for the transition line 
between phases O and II remains valid as r is increased, but for large enough r the transition 
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FIG. 7: The phase diagram in the K — 7 plane, with r = 2 > 1.75. Dashed and solid lines denote 
first- and second-order transitions respectively. The dotted line corresponds to the line 3TT 2 7 = 1. 




FIG. 8: Phase diagram in the -^00(7, r) plane. Solid lines represent continuous phase transitions, 
whilst the points indicate a first-order transition line. 
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FIG. 9: Density p plotted in the if 00(7,1") plane as a function of 7 for different values of r. 
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FIG. 10: Density p plotted in the K^ipj^r) plane as a function of r for different values of 7. 

is found to be first order, see the phase diagrams in figures [U] and [71 We may determine the 
condition for this change in behaviour by allowing the parameters x$ and X2 to be vanishingly 
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FIG. 11: Density p plotted in the K^ipj, r) plane as a function of t for different values of 7. 

small but non-zero along the transition line. We obtain the condition r = 3/2 . Thus the 
O-II line is second order and given by the equation ?>K 2 ^ = 1 for r < 3/2 and it is first 
order for r > 3/2, and falls below the line 3i^ 2 7 = 1. This behaviour may be understood in 
another way; as 7 — > 00, the walk becomes a doubled up self-avoiding walk, with an effective 
fugacity per bond K = K 2, y. We may then apply the known results to find that for K < 1/3 
the transition from phase O to the dense phase is continuous, and for K > 1/3 the transition 
is first-order. There is a transition when r = 3/2 and 3K = 3K 2/ y = 1. Since phase II 
is a saturated phase, the walk remains a doubled up self- avoiding walk for finite 7, hence 
the observed result, and the prediction that for r = 3/2 the transition line between phase 

and phase II is of the same type as the transition. It would be interesting to ascertain 
whether this prediction is maintained if a full calculation for the model is performed. For 
t < 3/2 the first-order line between phases O and I and the first-order line between phases 

1 and II are tangential to each other at the critical end point which defines the end of the 
critical line between phases O and II. However, for 3/2 < r < 1.75 (see figure [6]) the three 
first-order lines meet at a triple point (however the first-order lines between phases O and I 
and between I and II are still tangential). 

Another interesting change arises when r w 1.75 (see figure [7]). The phase transition 
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between phases I and II becomes totally second-order, and the first-order boundary to phase 
O becomes smooth, and the point where the three phases meet is again a critical end point. 




FIG. 12: Fraction of paired bonds plotted in the X 00 (7,r) plane as a function of r for different 
values of 7. 

Of particular interest is the phase diagram in the plane K^j, r), where the average 
length of the walk first diverges. This plane is particularly important, since it corresponds 
to the phase diagram which is seen in a canonical simulation of single finite-sized walk in 
the limit of infinite walk length. Indeed it is sometimes referred to as the "thermodynamic 
limit" for the canonical dilute polymer problem. The phase diagram in the plane is 
shown in figure [HI There are four different phases, differentiated by the values of p and 
$. There are two critical p = phases, both of which should be in the self-avoiding walk 
universality class. One has $ = (the standard SAW phase) whilst the other has $ = 1 
corresponding to a double SAW phase. There are two non-critical phases, one with $ = 1, 
corresponding to a standard type of collapsed phase as seen in the ©-point model, except 
that the walk is doubled. The other phase has a value of $ different from both and 1. 

It is interesting to note that 7e( / r) decreases as r increases, making it easier for the walk 
to pass from the p = 0, $ = phase to the neighbouring dense phase as r is increased. At 
first sight this seems odd, however the formation of some double bonds enables the formation 
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of a branched polymer type conformation, which is more naturally space filling. This would 
lead one to expect that perhaps the phase boundary between the p = 0, $ = phase and 
the p>0,0<$<l phase may not be in the G-point universality class. The mean-field 
nature of the Bethe approximation does not enable one to investigate this point. 

From previous discussion, the line separating the two $ = 1 phases, found at r = 3/2, 
may be expected to be a line of tricritical points in the 0-point universality class. This line 
appears to end at a tricritical end point on a first-order line (r = 3/2,7 w 5.89), which in 
its turn ends at another tricritical point at r rs 1-75,7 ~ 7.845. Whilst much could change 
once fluctuations are included, one could expect much of the behaviour seen here to remain 
in the real system. Figures [9] to [T2] show different plots of p and $ showing their behaviour 
at the different phase boundaries, confirming the orders of the transitions shown. 



V. DISCUSSION 



In this article we have investigated the phase diagram of a non-crossing random walk 
model on the square lattice where the walk is allowed to visit lattice bonds twice, and 
the site configurations are chosen to represent the zipping/unzipping of a DNA molecule. 
Solvent quality is included through attractive nearest-neighbour interactions. This model 
is examined in the framework of the Bethe approximation and found to have a rich phase 
diagram. Whilst there are limitations to the method, it usually captures the essential features 
of the phase diagramjjj], [ljj. It would be of interest to look at this model using some other 
method to confirm and examine further the phase diagram presented. 
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